
#include <boost/date_time/posix_time/posix_time.hpp>
#include <boost/filesystem.hpp>
#include <boost/format.hpp>
#include <stdlib.h>
#include <time.h>
#include "ADefine.h"
#include "ATimeSpace.h"
#include "ACatUCAC4.h"
#include "SimStar.h"
#include "OrbitDef.h"

using namespace boost::posix_time;
using namespace boost::filesystem;
using namespace AstroUtil;

// 在UCAC4星表中依据搜索条件查找最亮星
// (ra0, dec0): 中心坐标, 角度
// radius: 半径, 角分
// (ra, dec): 亮星J2000坐标
bool FindStar(double ra0, double dec0, double radius, double &ra, double& dec, double &mag0) {
	ACatUCAC4 ucac4;
    int i, n;

	ucac4.SetPathRoot("/Users/lxm/Catalogue/UCAC4");
	if (ucac4.FindStar(ra0, dec0, radius)) {
		ucac4item_ptr stars, star, t = nullptr;
        short mag(30000);

		stars = ucac4.GetResult(n);
		for (i = 0, star = stars; i < n; ++i, ++star) {
            if (star->apasm[1] > 5000) continue;
            if (star->apasm[1] < mag) {
                mag = star->apasm[1];
                t = star;
            }
		}
        if (t == nullptr) return false;
        ra = double(t->ra) / AU_MAS;
        dec= double(t->spd) / AU_MAS - 90.0;
        mag0 = t->apasm[1] * 0.001;

        return true;
	}
    return false;
}

void UTC2SatTime(ptime& utc, OrbitNode& node) {
    ptime bjt = utc + hours(8);
    ptime bjt0 = from_iso_extended_string("2024-09-21T18:00:00");
    int st0 = 54381600;
    ptime::date_type today = bjt.date();
    ptime::time_duration_type hms = bjt.time_of_day();
    int YY = today.year();
    int MM = today.month().as_number();
    int DD = today.day();
    int hh = hms.hours();
    int mm = hms.minutes();
    int ss = hms.seconds();

    node.tmSat = st0 + (bjt - bjt0).total_seconds();
    snprintf (node.bjtstr, 40, "%d%02d%02d-%02d:%02d:%02d.000",
        YY, MM, DD, hh, mm, ss);
}

// 生成时段内仿真轨道
void GenerateSimStar(const char* bjtBeg, const char* bjtEnd, int period) {
    srandom((unsigned int) time(NULL));

    ptime bjt0 = from_iso_extended_string(bjtBeg);
    ptime bjt1 = from_iso_extended_string(bjtEnd);
    int periodSec = period * 60;
    // int periodSec = 7200; // Vega...
    int YY, MM, DD, hh, mm, ss;
    int code(0);
    int date = 0;

    double azi0 = 180 * AU_D2R;
    double alt0 = 60 * AU_D2R;
    double azi, alt;
    double ha, raApp, decApp;
    double ra, dec, mag;
    double radius(10 * 60.0); // 5度;半径
    ATimeSpace ats;
    ats.SetSite(116.26663, 40.06856, 88, -8);

    OrbPtrVec orbs;
    // 生成轨道
    printf ("generate orbit for star....\n");
   for (ptime bjt = bjt0; bjt < bjt1; bjt += minutes(period)) {
        azi0 = (random() * 3200.0 / RAND_MAX * 0.1 + 20.0) * AU_D2R;
        alt0 = (random() * 450.0 / RAND_MAX * 0.1 + 35.0) * AU_D2R;

        // 生成中心赤道坐标
        ptime utc0 = bjt - hours(8);
        ats.SetMJD(utc0.date().modjulian_day()
            + utc0.time_of_day().total_seconds() / 86400.0);
        ats.Horizon2Eq(azi0, alt0, ha, decApp);
        raApp = cycmod(ats.LocalSiderealTime() - ha, AU_2PI);
        if (!FindStar(raApp * AU_R2D, decApp * AU_R2D, radius, ra, dec, mag)) {
            printf ("skipped !!!\n");
            continue;
        }
        printf ("\t%s: AA=%8.4lf %8.4lf; EQ=%8.4lf %8.4lf; Mag=%.3lf\n",
            to_iso_extended_string(bjt).c_str(),
            azi0 * AU_R2D, alt0 * AU_R2D,
            ra, dec,
            mag);
        // Vega...
        // ptime utc0 = bjt0 - hours(8);
        // ptime bjt = bjt0;
        // Vega...
        ptime::date_type today = bjt.date();
        ptime::time_duration_type hms = bjt.time_of_day();
        YY = today.year();
        MM = today.month().as_number();
        DD = today.day();
        hh = hms.hours();
        mm = hms.minutes();
        ss = hms.seconds();

        OrbitPtr orbit(new Orbit);
        orbit->date = YY * 10000 + MM * 100 + DD;
        orbit->time = hh * 10000 + mm * 100 + ss;
        if (date != orbit->date) {
            code = 0;
            date = orbit->date;
        }
        orbit->code = ++code;

        OrbNodeVec &nodes = orbit->nodes;
        // Vega
        // ra = (18 + (36 + 56.34 / 60) / 60.0) * 15.0;
        // dec= 38 + (47 + 1.28 / 60) / 60.0;
        // Vega
        for (int i = 0; i < periodSec; i += 5) {
            OrbitNode node;
            ptime utc = utc0 + seconds(i);
            ats.SetMJD(utc.date().modjulian_day() + utc.time_of_day().total_seconds() / 86400.0);
            ats.EqTransfer(ra * AU_D2R, dec * AU_D2R, raApp, decApp);
            ats.Eq2Horizon(ats.LocalSiderealTime() - raApp, decApp, azi, alt);

            UTC2SatTime(utc, node);
            node.ra  = ra;
            node.dec = dec;
            node.azi = azi * AU_R2D;
            node.alt = alt * AU_R2D;
            nodes.push_back(node);
        }
        orbs.push_back(orbit);
    }

    // 存储轨道
    printf ("output >> single orbit >> begin\n");
    path pathDir;
    date = 0;
    for (auto &it : orbs) {
        if (it->repeated) continue;

        if (it->date != date) {
            date = it->date;
            pathDir = path((boost::format("STAR%d") % date).str());
            if (!exists(pathDir)) create_directories(pathDir);
        }
        boost::format fmt("%dT%06d_%d.txt");
        fmt % it->date % it->time % it->code;
        path pathName(pathDir);
        pathName /= fmt.str();
        printf ("\t %s\n", fmt.str().c_str());

        FILE *fp = fopen(pathName.c_str(), "w");
        it->Output(fp);
        fclose(fp);
    }
    printf ("output >> single orbit >> over\n");

    printf ("Over\n");
}
